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We discuss simulations of finite temperature nuclear matter on the lattice. We introduce a 
new approximation to nucleon matrix determinants that is physically motivated by chiral effective 



OO , 

I theory. The method involves breaking the lattice into spatial zones and expanding the determinant 

CN ' in powers of the boundary hopping parameter. 
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I. INTRODUCTION 



We consider quantum simulations of nuclear matter on the lattice. In particular we ad- 
dress the problem of calculating the contribution of nucleon/nucleon-hole loops at nonzero 
nucleon density. With the help of auxiliary boson fields, all nucleon interactions can be 
written in terms of one-body interactions in a fluctuating background. In the grand canoni- 
cal ensemble, the contribution of nucleon/nucleon-hole loops to the partition function equals 
the determinant of the one-body interaction matrix. Since the determinant of the interac- 



tion matrix for a general 



boson field configuration is not positive, stochastic methods such 



as hybrid Monte Carlo do not give the sign or phase of the determinant. Instead one 

must rely on much slower and more memory intensive algorithms based on LU factorization, 
which decomposes matrices in terms of a product of upper and lower triangular matrices. 

The number of required operations in LU factorization for an n x n matrix scales as n^. 
It has been shown in the literature that repeated calculations of matrix determinants with 
o„l. ■ocaU.ed c.a,,. can be ..a,.li„ed in vanon. wa.= M- Howeve. it i. difficult to 
avoid the poor scaling inherent in the method. If V is the spatial volume and (3 is the inverse 
temperature measured in lattice units, a simulation that includes nucleon/nucleon-hole loops 
requires times more operations than the corresponding quenched simulation without 

loops. This slowdown should not be confused with the infamous fermion sign or phase 
problem which becomes significant at temperatures T < 1 MeV. The computational 
bottleneck we are considering is due to the inefficiencies of the algorithm and persists at 
all temperatures. It is this numerical challenge which sets current limits on nuclear lattice 
simulations. 

In this paper we introduce a new approach to approximating nucleon matrix determinants. 
We begin with a review of the current status of nuclear matter simulations on the lattice and 
look to chiral effective theory to determine the relative importance of various interactions. 
We then introduce the concept of spatial zones and suggest a new expansion of the nucleon 
determinant in powers of the hopping parameter connecting neighboring zones. Rigorous 
bounds on the convergence of this expansion are given as well as an estimate of the required 
size of the spatial zones as a function of temperature. We apply the expansion to a realistic 
lattice simulation of the interactions of neutrons and neutral pions. 
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II. NUCLEAR LATTICE SIMULATIONS 

Recently Miiller, Koonin, Seki, and van Kolck 0], pioneered the study of quantum many 
body effects in infinite nuclear matter at finite density and temperature. In their work they 
considered only nucleon degrees of freedom and used an effective Hamiltonian of the form 

H = K+V, + V„, (1) 

where K is the kinetic energy, Vc is the two-body scalar potential, and is the two- 
body tensor potential. While the simulation was done on the lattice we will write their 
Hamiltonian in the more familiar continuum language. In the continuum Vc-, and V^r 



take the form 

K = -^ j d'x^lV'^^^, (2) 



In our notation summations are implied over repeated indices. is the nucleon mass. 

ipl^lipar) creates(annihilates) a nucleon of spin a and isospin r, and ct^tkA are the elements 
of a generalized Pauli spin-isospin matrix. Both potentials are assumed to have Skyrme-like 
on-site and next-nearest-neighbor interactions. 



Vcix - x') = V^^^Six - x') + V}^^V^5ix - f ), (5) 

V,{x - x') = V}^^6{x - x') + Vj^^V^6{x - x'). (6) 

The grand canonical partition function is given by 

Z = Tr [exp [-(3 {H - /i.n.)]] , (7) 



where /Zt- is the isospin-dependent chemical potential and n,- is the nucleon number operator 
for isospin index r. We can rewrite the quartic interactions in Vc and V^ using the Hubbard- 
Stratonovich transformation The Hubbard-Stratonovich transformation uses the 

3 



identity, 
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where A is any quantum operator. This allows the mapping of the interacting nucleon 
problem to a system of noninteracting nucleons coupled to a fluctuating background field. 
With this transformation the expectation value of any observable O can be written as 

^ /D0G(0)det(M(0))O(0) 
^ ' J D<pG{<f))det{M{(f))) ' ^ 

where collectively represents the Hubbard-Stratonovich fields (as well as any other bosonic 
fields), M(0) is the one-body nucleon interaction matrix, and G(0) is a function of the 0's. 

Using this formalism Miiller, et. ai, were able to measure the thermodynamic properties 
of nuclear matter and find signs of a first order phase transition from an uncorrelated Fermi 
gas to a clustered phase. They examined temperatures from 3.0 MeV to 100 MeV with up to 
30 time steps and a spatial volume of 4'^ with lattice spacing a = 1.842 fm. Unfortunately 
the LU factorization algorithm used to compute determinants in the simulation scales as 
[VPY and thus going beyond lattice systems of size 4^ is problematic. 

III. CHIRAL EFFECTIVE THEORY AND NUCLEAR FORCES 

There have been several recent efforts to describe nuclear forces starting from chiral 
effective theory. This line of study was initiated by Weinberg 3Q[ll|. The idea is to 
expand the nuclear interactions in powers of q/A^, where q is the typical external momentum 
of the nucleons and is the chiral symmetry breaking scale or equivalently the hadronic 
mass scale (~ 1 GeV). The momentum cutoff scale A for the effective theory is set below 
A^, and the renormalization group flow of operator coefficients from scale A^ to A suppress 
the effects of higher dimensional operators. Chiral symmetry in the limit of zero quark 
mass imposes additional constraints on the possible momentum and spin dependence of the 
interaction terms. Assuming naturalness for the renormalized coupling constants in the 
Lagrangian at the scale A^, one expects in the low energy effective theory that contributions 
to nucleon forces from operators at higher order in the chiral expansion are negligible. 

Weinberg's work was followed by applications of chiral effective theory to the nucleon 



potential 



12| and alternative approaches to power counting without apparent fine tunin g in 



the presence of long scattering lengths 



3 to power countmg witnout apparent nne tunmgi 
QQ. Recent low energy studies QQHE 



have also integrated out pion fields to produce energy independent two- and three- nucleon 
potentials, and the effective theory approach has been used to calculate nuclear spectra as 
well as phase shifts and scattering lengths which compare favorably with potential model 
calculations. 

In Weinberg's power counting scheme one deals with infrared singularities in bound state 
problems by distinguishing between reducible and irreducible diagrams. Reducible diagrams 
are those that can be disconnected by cutting internal lines that correspond with particles 
in the initial or final state. In the notation of the power of Q'/A^ for any irreducible or 
non-reducible diagram is given by 



V 



A-^ + 2L-2C + Y,yA 



(10) 



where Ef is the number of external nucleon lines, L is the number of loops, C is the number 
of connected pieces, Vi is the number of vertices of type and 5i is the index of vertex i. 
The index 5i is given by 

5^ = cii + f - 2, (11) 

where di is the number of derivatives and fi is the number of nucleon fields in the vertex. 
We let N represent the nucleon fields. 



N 



p 




r 




n 







(12) 



We use Tj to represent Pauli matrices acting in isospin space, and we use a to represent Pauli 
matrices acting in spin space. Pion fields are notated as vTj, and /x is the nucleon chemical 
potential. We denote the pion decay constant as ^ 183 MeV and let 



D = l + r^j/Fl 



(13) 



The lowest order Lagrange density for low energy pions and nucleons is given by terms with 
5i = 0, 



£(0) 



-1 r-l 



(Vtt, 



^2 ^2 

— TT, 



Tiff ■ VvTi 



- iD'^mln^ + N[ido - (m^ - i2)]N 
N - D-'F-^N[e,,kr^7T^7^k]N 



- ICsNNNN - ^CrNaN ■ NaN. 



(14) 
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gA is the nucleoli axial coupling constant, and eijk is the Levi-Civita symbol. For the purposes 
of power counting, the pion mass m-Tr is equivalent to one power of the momentum q. The 
nucleon mass term NN actually has index 6i = —1. However the coefficient of this term is 
fine-tuned using fi to set the nucleon density, and so the NN term is reduced to the same 
size as other terms with index 6i = 0. At next order we have terms with 6i = 1, 

^'•'^ = 2^NV'N + ... (15) 

The important point for our discussion is that the lowest order Lagrange density C^'^^ de- 
scribes static nucleons. Spatial hopping of nucleons first appears at subleading index 6i = 1. 
This suggests that in some cases one could compute the determinant of the nucleon inter- 
action matrix as an expansion in powers of the spatial hopping parameter. We should note 
one point of caution though. The C^^^ kinetic energy term cannot be ignored if the infrared 
singularities of C^^^ are not properly dealt with. Since diagrammatic methods are not used 
in nuclear matter Monte Carlo simulations, we cannot separate reducible and irreducible 
diagrams. However if the simulation is done at non-zero temperature T then that will serve 
to regulate the infrared singularities. We will explicitly see the effect of temperature on the 
convergence of the expansion later in our discussion. 



IV. SPATIAL ZONES 

In spatial zones were used to calculate the chiral condensate for massless QED in 
three dimensions. In that paper the main problem was dealing with sign and phase problems 
that arise in the Hamiltonian worldline formalism with explicit fermions. The idea of the 
zone method is that fermions at inverse temperature (3 with spatial hopping parameter h 
have a localization length of 

/ ~ V^- (16) 

We now apply the zone idea to our determinant problem. Let M be the nucleon matrix, 
in general an n x n complex matrix. We partition the lattice spatially into separate zones 
such that the length of each zone is much larger than the localization length /. Since most 
nucleon worldlines do not cross the zone boundaries, they would not be affected if we set 
the zone boundary hopping terms to zero. Hence we anticipate that the determinant of M 
can be approximated by the product of the submatrix determinants for each spatial zone. 
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Let us partition the lattice into spatial zones labelled by index j. Let {Pj} be a complete 
set of matrix projection operators that project onto the lattice sites within spatial zone j. 
We can write 

M = J2 PjMP, = Mo + Me, (17) 

id 

where 

Mo = 5^P.MP„ (18) 

i 

If the zones can be sorted into even and odd sets so that 

PjMPi = (20) 

whenever i is even and j is odd or vice- versa, then we say that the zone partitioning is 
bipartite. We now have 

det(M) = det(Mo) det(l + M^^Me) 

= det(Mo) exp (trace (log (1 + Mq^Me))). (21) 

Using an expansion for the logarithm, we have 

det(M) = det(Mo) exp tTeice{{M^^MEf)j • (22) 

Let us define 

A™ = det(Mo) exp V ^ trace((Mo-^Mij)P) . (23) 

\U P J 

Let Xk{MQ^ Me) be the eigenvalues of Mq^Me and R be the spectral radius, 

R= max {\XkiM^^ME)\). (24) 

k=l,...,n 

It has been shown [2^ that for R < 1 

|det(M)-A^| ^^^^ 



< cP^^e'^" (25) 



where 

c= -nlog(l-i?). (26) 
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The spectral radius R determines the convergence of our expansion. R can be reduced by 
increasing the size of the spatial zone relative to the localization length I. In the special 
case where the zone partitioning is bipartite, we note that for any odd p, 

trace((Mo-^Ms)P) = 0. (27) 

In that case 

^2m+l = ^2m, (28) 

and so we gain an extra order of accuracy without additional work. 

V. APPLICATION TO NEUTRON MATTER SIMULATIONS 

We now illustrate the zone determinant expansion for a simple but realistic lattice sim- 
ulation of neutron matter. The formalism we use is a merger of chiral effective theory and 
Euclidean lattice methods. In our analysis we focus on the convergence and accuracy of 
the zone determinant expansion method. In order to provide a detailed quantitative as- 
sessment we compare the zone determinant approximations with exact determinant results. 
Given the severe limitations of the exact determinant method, we are not able to probe 
large volumes nor comment on finite volume scaling. However a full account of results for 
large volume simulations as well as an introduction to the new approach combining chiral 



effective theory and lattice methods is forthcoming in another paper [21 1. 

Our starting point is the same as that of Weinberg 9]. We start with the most general 
local Lagrange density involving pions and low-energy nucleons consistent with Lorentz 
and translational invariance, isospin symmetry, and spontaneously broken chiral symmetry. 
This will produce an infinite set of interaction terms with increasing numbers of derivatives 
and/or nucleon fields. The dependence of each term on the pion field is governed by 
the rules of spontaneously broken chiral symmetry. Degrees of freedom associated with 
heavier mesons such as the p, heavier baryons such as the A's as well as antinucleons are all 
integrated out. We also integrate out nucleons with momenta greater than the cutoff scale 
A. The contribution of these effects appear as coefficients of local terms in our pion-nucleon 
Lagrangian. 

For simplicity we consider only neutrons and neutral pions. We let ip represent the 
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neutron spin states, 



The terms in our effective pion-nucleon Euclidean action are 



(29) 



where 



NNNNj 



(30) 



(31) 



Snn = Jd'rdr, [^^^ - + (m^ - /x)V^t^ 



(32) 



■NN 



-gijj^ ( a ■ Vtto ) 



(33) 



Snnnn = I d^rdr^ [|C:V^V^V:] • (34) 

We have kept terms in the lowest order chiral Lagrange density C^^^ containing neutrons 
and neutral pions. We have dropped the factors of D^^ which at lowest order contribute to 
multipion processes. We have also chosen to include the neutron kinetic energy term even 
though it appears in C^^\ This is useful if we wish to recover the exact free neutron Fermi 
gas in the weak coupling limit. 

On the lattice we let a be the spatial lattice spacing and at be the temporal lattice 
spacing. We use the notation 1,2,3, or 4 to represent vectors that extend exactly one 
lattice unit (either a or a^) in the respective direction. We use dimensionless lattice fields 
and dimensionless masses and couplings by multiplying by the corresponding power of a. 
For example, rh.„ = ■ a, rhiy = niiy ■ a, g = g ■ a^^, ft = fi ■ a, and C = C ■ a^^. We use 



standard fermion path integral conventions at finite time step 

ip'{n) = ilj{n — 4) 



22 



and define 



(35) 



in order to write the lattice path integral in standard form. The lattice actions have the 
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form 



n n,/=l,2,3 



(36) 



ft n,/=l,2,3 



in 



(37) 



^-^iv = - + 3) - vr(n - 3)]] 

n 

- [vr(n + i) -7r(n- i) -z7r(fi + 2) +i7r(n-2)]] 

n 

- [7r(n + i) -7r(n- i) +i7r(n + 2) -i7r(n-2)]] 



(38) 



S 



NNNN 



(39) 



We can reexpress Sfj]\[f^]\[ using a discrete Hubbard-Stratonovich transformation ^] for 
(7 < 0, 



exp(-^^^|(n)^;(n)^t(n)V^i(n)) 



i 5^ exp \-{^ + Xs){^^n)^\.{n) + V^I(n)V^i(n) - 1) 



s=±l 



where A is defined by 



cosh A = exp(— — 



2a 



(40) 



(41) 



VI. RESULTS 



For the simulation results presented in this section we use the values a =150 MeV, 
a^^ = 225 MeV, C = -4.0 ■ 10^^ MeV'^, and ^ = 6.8 • 1^^ MeV^^ The value of g is set 
according to the tree level Goldberger-Treiman relation 



g = P-^QA = 6.8 ■ 10-3 MeV"\ 



(42) 
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The value of C is tuned to match the s-wave phase shifts on the lattice for nucleon scattering 
at lattice spacing (150 MeV)^^. The calculation of phase shifts on the lattice is discussed 
in a forthcoming article which details the entire formalism j2l| . 

We present data for three independent pion and Hubbard-Stratonovich field configura- 
tions on a 4^ X 6 lattice at temperature T = 37.5 MeV and neutron density p = 0.57pnuci 
where nuclear density is 

Pnuci = 2.8 ■ 10^^ g/cin = 1.2 ■ 10^ MeV^ (43) 

We refer to spatial zones according to their x,y,z lattice dimensions [n^;, n^, n^]. At this 
lattice spacing and temperature our localization length estimate is 

/ ^ = 0.57, (44) 

and so we expect the zone determinant expansion to converge for even the smallest zones, 
[nx,ny,nz] = [1,1,1], consisting of groups of lattice points with the same spatial coordi- 
nate. In Table 1 we show the spectral radius R of Mq^Me and the determinant expansion 
for [1, 1, 1] zones for three independent pion and Hubbard-Stratonovich configurations at 
equilibrium. 



configuration 


#1 


#2 


#3 


R 


0.538 


0.523 


0.534 


log(det(M)) 


16.4727 + 0.3666i 


18.1193 + 0.4479i 


18.2612 -0.1811i 


log(Ao) 


12.1700 + 0.3807i 


13.7933 + 0.4900Z 


14.0060 - 0.2075i 


log(A2) 


16.4964 + 0.3686i 


18.1792 + 0.44772 


18.2841 -0.1801i 


log(A4) 


16.4959 + 0.3659i 


18.1329 + 0.4468i 


18.2856 -0.1805i 


log(A6) 


16.4664 + 0.3667i 


18.1147 + 0.4482Z 


18.2548 -0.1813i 


log(A8) 


16.4739 + 0.3666i 


18.1203 + 0.4478i 


18.2622 -0.1810i 



Table 1: Convergence of determinant series for [1, 1, 1] zones 



We see that the spectral radius R is less than 1 as expected from the localization length 
estimate. We also find that the rigorous bounds in (j25j) are satisfied. Since most of the 
eigenvalues of Mq^Me are much smaller in magnitude than the spectral radius R, we observe 
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that the prefactor c in ()25|) is actually much larger than needed for these examples. From 
the data in Table 1 we find empirically 



|det(M)-A2„| |det(M)-A 



2m+l I 



R 



2m+l 



lA I lA I " ■ 

\^2m\ |^2m+l| 

We now repeat the same analysis with the same lattice dimensions and temperature 
but at much higher density, p = 1.67pnuci- The results for three independent pion and 
Hubbard-Stratonovich field configurations are shown in Table 2. 



configuration 


#1 


#2 


#3 


R 


0.520 


0.478 


0.536 


log(det(M)) 


70.6349 + 0.4121i 


76.1762 + 0.7098i 


73.4831 -0.1576i 


log(Ao) 


64.8000 + 0.4406i 


71.0009 + 0.7989i 


67.8503 - 0.1809i 


log(A2) 


70.9530 + 0.4120i 


76.3927 + 0.7005i 


73.7923 -0.1538i 


log(A4) 


70.6032 + 0.4114i 


76.1605 + 0.7106Z 


73.4483 - 0.1583i 


log(A6) 


70.6390 + 0.4124i 


76.1778 + 0.7097i 


73.4888 -0.1574i 


log(A8) 


70.6343 + 0.4120i 


76.1760 + 0.7098Z 


73.4819 - 0.1576i 



Table 2: Determinant series for [1,1,1] zones at p = 1.67pnuci 
Comparing Tables 1 and 2, we conclude that the zone determinant expansion appears to be 

unaffected by the increase in nucleon density. 

On the other hand as the temperature decreases, the localization length increases and 

therefore the convergence of the zone determinant expansion should slow down. In Table 3 

we show the expansion for a 6'^ x 6 lattice at T = 37.5 MeV, 6^ x 9 lattice at T = 25.0 MeV, 

6^ X 12 lattice at T = 18.8 MeV, 6^ x 18 lattice at T = 12.5 MeV. The chemical potential 

is kept at the same value used for the results in Table 1 (p = 0.8mAr). 



T 


37.5 MeV 


25.0 MeV 


18.8 MeV 


12.5 MeV 


R 


0.5122 


0.6742 


0.7631 


0.8638 


log(det(M)) 


65.8009 - 0.7250i 


37.6912 - 1.5930i 


16.0023 + 0.0541i 


5.7618 + 0.0452i 


log(Ao) 


51.3988 - 0.7888i 


20.7653 - 1.57522 


4.1999 + 0.2033i 


0.1796 - 0.0019i 


log(A2) 


66.0028 - 0.7218i 


36.2570 - 1.6535Z 


11.7590 + 0.0861i 


1.5071 + 0.0222i 


log(A4) 


65.8289 - 0.7243i 


38.0872 - 1.5800i 


15.9390 + 0.0337i 


3.8579 + 0.0399i 


log(A6) 


65.7926 - 0.7253i 


37.6528 - 1.5924i 


16.3121 + 0.0478i 


5.5534 + 0.0475i 


log(A8) 


65.8026 - 0.7249i 


37.6780 - 1.5943i 


16.0111 + 0.0575i 


6.0067 + 0.0473i 
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Table 3: Determinant series for [1, 1, 1] zones for varying temperatures 



The increase in R and the slowdown in the series convergence is consistent with our intuition 
based on the localization length. Although we are looking at neutron matter in this example, 
we comment that this is the appropriate temperature range for observing the liquid-gas 
transition in symmetric nuclear matter (T ~ 20 MeV) [2^[23[3|- The localization length 
becomes greater than 1 for T ~ 10 MeV, and so we expect the series to break down for 
[1, 1, 1] zones at colder temperatures. This also happens to be the temperature at which long 
range interactions due to pairing become important. In order to see pairing correlations, 
one should therefore use larger zone sizes. However we note that at temperatures much less 
than 10 MeV there is also a significant sign problem, and the numerical difficulties will be 
substantial independent of the method used to calculate determinants. 

We now look at how convergence of the expansion is improved by using larger spatial 
zones. In Table 4 we show the determinant expansion for a 6^ X 6 lattice at T = 37.5 MeV 
for [1,1,1], [2,2,2], and [3,3,3] zones H- 



zone 


[1,1,1] 


[2,2,2] 


[3,3,3] 


R 


0.5122 


0.3013 


0.3014 


log(det(M)) 


65.8009 - 0.7250i 


65.8009 - 0.7250i 


65.8009 - 0.7250i 


log(Ao) 


51.3988 - 0.7888i 


58.6585 - 0.75432 


61.0471 - 0.7722i 


log(A2) 


66.0028 - 0.7218i 


65.8570 - 0.7231Z 


65.8317 - 0.7233Z 


log(A4) 


65.8289 - 0.7243i 


65.8015 -0.7251i 


65.8007 - 0.7250i 


log(A6) 


65.7926 - 0.7253i 


65.8009 - 0.72502 


65.8009 - 0.7250i 


log(A8) 


65.8026 - 0.7249i 


65.8009 - 0.7250i 


65.8009 - 0.72502 



Table 4: Determinant series for various zone sizes 

It is somewhat unusual that the spectral radius R is about the same for [2, 2, 2] and [3, 3, 3], 
however the convergence of the series clearly improves as we increase the zone size. 

We now investigate the convergence of the determinant zone expansion for physical ob- 
servables. Let us return to the data in Table 1 for a moment. We observe that at any given 
order the error appears to be about the same for each of the three independent configura- 
tions. Since the measurement of a physical observable does not dependent on the overall 
normalization of the partition function, this suggests that the zone determinant expansion 
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<p I (r)p^(0)> exact determinant 




3 3 z 

y 



FIG. 1: Opposite spin radial distribution function using exact determinants. 

could be more accurate in approximating physical observables. We now check to see if this 
is so for a particular example. 

Let us define the neutron occupation number at site f, 

PTa)(0 = ^Ta)(O^Ta)(0 = ^Ta)(^)^Ta)(^+ 4)- (46) 

We measure the opposite spin radial distribution function, 

<Pt(r)Pi(0)>, (47) 

by sampling 20 independent pion and Hubbard- Stratonovich field configurations. For = 
we plot the results for the opposite spin radial distribution function in Fig. 1 using exact 
matrix determinants for a 4'^ x 6 lattice at T = 37.5 MeV. In Figs. 2-6 we show the 
error in the radial distribution function if the estimate is used in place of det(M) for 
m = 0,2,4,6,8. 

The approximation at order is a lot better than expected given the error of Aq in ap- 
proximating det(M). Overall we find that the zone expansion is significantly more accurate 
for the radial distribution function than the expansion of the determinants. It is premature 
to say if this is typical of all physical observable measurements. Nevertheless if most of the 
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error of is in fact independent of the pion and Hubbard-Stratonovich configurations at 
equilibrium, tlien we expect an improvement in accuracy for most physical observables. 

VII. SUMMARY AND CONCLUSIONS 

We have discussed lattice simulations of finite temperature nuclear matter and a new 
approximation method called the zone determinant expansion for nucleon matrix determi- 
nants. The expansion is made possible by the small size of the spatial hopping parameter. 
We know from power counting in chiral effective theory that the spatial hopping parameter is 
suppressed relative to the leading order interactions at low energies. The zone determinant 
expansion is given by 

det(M) = det(Mo) exp tmce{{MQ^MEy) j , (48) 

where M is the nucleon one-body interaction matrix, Me is the submatrix consisting of zone 
boundary hopping terms, and Mq is the submatrix without boundary hopping terms. The 
convergence of the expansion is controlled by the spectral radius of Mq^Me- Physically we 
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Error in <p , (r)p^(0)>, order=2 



X 10 



1.5-, 




0.5- 
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y 

FIG. 3: Error in the opposite spin radial distribution function at order m = 2. 
expect the convergence to be rapid if the locahzation length of the nucleons 



is small compared to the size of the spatial zones. 

We tested the zone determinant expansion using lattice simulations of neutron matter 
with self interactions and neutral pion exchange. The convergence of the expansion was 
measured for several configurations at temperature T = 37.5 MeV and using [1, 1, 1] spatial 
zones. By decreasing the temperature from T = 37.5 MeV to 12.5 MeV we found that the 
convergence of the expansion becomes slower, as predicted by the increase in the localization 
length. But we then showed that convergence could be accelerated by increasing the size 
of the zones from [1, 1, 1] to [3, 3, 3]. Finally we looked at the convergence of the expansion 
for the opposite spin radial distribution function 



We found that the accuracy of the expansion for this physical observable was significantly 
better at each order than that for the expansion of the determinants. 

The number of required operations for calculating the nucleon determinant using LU 
factorization for an n x n matrix scales as n^. Therefore a nuclear lattice simulation that 




(49) 



<PT(^)Pi(0) > • 



(50) 
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includes nucleon/nucleon-hole loops requires {VPY times more operations than the quenched 
simulation without loops. This numerical challenge has been the most pressing limitation 
on finite temperature nuclear lattice simulations to date. 

For the zone determinant expansion method at fixed zone size, the computation cost 
scales only as f{m)(3^ where m is the order of the expansion. For a simulation on a lattice 
with spatial dimensions 8^, one can accelerate the simulation by a factor of about 10^ to 
10^, depending on the expansion order and size of the spatial zone. The savings are greater 
on larger lattices and should facilitate future work in the area of finite temperature nuclear 
lattice simulations. 
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